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Abstract. The spatially localized bound states of two electrons in the adiabatic two-dimensional Holstein- 
Hubbard model on a square lattice are investigated both numerically and analytically. The interplay 
between the electron-phonon coupling g, which tends to form bipolarons and the repulsive Hubbard inter- 
action V > 0, which tends to break them, generates many different ground-states. There are four domains 
in the g, v phase diagram delimited by first order transition lines. Except for the domain at weak electron- 
phonon coupling (small g) where the electrons remain free, the electrons form bipolarons which can 1) 
be mostly located on a single site (small v, large g); 2) be an anisotropic pair of polarons lying on two 
neighboring sites in the magnetic singlet state (large v, large g); or 3) be a " quadrisinglet state" which 
is the superposition of 4 electronic singlets with a common central site. This quadrisinglet bipolaron is 
the most stable in a small central domain in between the three other phases. The pinning modes and the 
Peierls-Nabarro barrier of each of these bipolarons are calculated and the barrier is found to be strongly 
depressed in the region of stability of the quadrisinglet bipolaron. 

PACS. 71.10.Fd Lattice fermion models (Hubbard model, etc.) - 71.38.-|-i Polarons and electron phonon 
interactions - 74.20.Mn Nonconventional mechanisms (spin fluctuations, polarons and bipolarons, res- 
onating valence bond model, anyon mechanism, marginal Fermi liquid, Luttinger liquid, etc.) - 74.25.Jb 
Electronic structure 
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to a quantum field of phonons. It has been well-known for quite unreasonable to expect the bipolarons to become 

several decades that when the electron-phonon couphng superfluid at a non-negligible temperature. This aspect 

increases too much, the BCS theory breaks down because of the problem has been emphasized recently in ref . pl| . 

of lattice instabilities ||^. As a consequence, rather low However, the argument used by these authors was based 

critical temperatures 30/^) were predicted as the upper on standard considerations that did not take into account 

bound for real BCS superconductors |3|. Many theories the effect of mass reduction we shall discuss in this and a 

have subsequently been developed to describe the strong subsequent paper p7| . 

couphng regime with the hope to predict the existence 

„ T-.^r. 1 . -,1 1 • 1 •,• 1 ; Indeed, in realistic physical models, the characteristic 

of non-BCS superconductors with high critical tempera- 

,, ,. , „, ifinni energy of the bare electrons is usually a few eV and is 

ture. After the discovery by Bednorz and MuUer |^,||,|6[ ^ ^ 

. • 1 1 • 1 1 1 • much larger than the phonon energies which is at most 

of cuprate materials, which can be superconducting at 

... -, 11-1 about a tenth of an eV. As a result, the quantum fluctu- 

temperatures as nigh as lUUiv or more, the bipolaron ap- 

, , , s . , , . , , Oi ations of the phonons become generally negligible as soon 

proach (among others) regained much interest |71|. 

as the electron-phonon coupling is strong enough to gen- 

Since Landau H , it has been acknowledged that a sin- ^ ,• i m ^ i • ^ ^- i ^ 

" erate bipolarons. I hen the potential interactions between 

gle electron (or cquivalently a pair of noninteracting elec- , , , . , , , , , , , . , , • ; • 

V i ./ i c3 ^j-^g bipolarons are much larger than their quantum kinetic 

trons coupled to a deformable classical field) may localize , . • ; , i • i ; ; 

energy, m that situation, the many bipoiaron structures 

in the potential created self-consistently by a deformation , , , , , •, , , rr • t • i • 

should be well described by an effective Ismg pseudospm 

of the field. The resulting object is called "polaron" for one tt -i. • i- .• ■ i .• t->- i ^i 

Hamiltoman, predicting an insulating Bipoiaron Charge 

electron or "bipoiaron" for two electrons. The bipoiaron 



Density Wave at low temperature |l2| , |l3yi4[| . 
theory of Alexandrov et al.|^ involves small bipolarons 

which are pairs of electrons with opposite spins, sharply However, there might exist special and exceptional sit- 

localized at single sites of the lattice. Actually, because nations where the effective mass of the bipolarons is not 

the phonons are quantum, these bipolarons are hard-core quasi-infinite but becomes small enough so that they pos- 

bosons that could condense in a superfluid state. For mod- sibly condense into a superfluid state. The smaller the 

els in two dimensions and more, bipolarons exist only bipoiaron mass is, the higher the critical temperature should 

when the electron-phonon coupling is large enough ]To| , be. As conjectured in ref.[p^ and [ p^ , this situation might 

and they are always sharply localized as small bipolarons be produced by a well-balanced interplay between the 

when the interactions are local. Thus, taking physically re- bare electronic kinetic energy, the electron-phonon cou- 

alistic parameters for the model, the effective mass of the pling and the direct electron-electron repulsion. The aim 

bipolarons becomes so huge (quasi-infinite) that it seems of this paper is to study this interplay in the simplest 
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Holstein-Hubbard (HH) model where these interactions The behavior of the bipolaron in the two-dimensional 
are present. case is quite different from the one-dimensional case. Al- 

though it does not describe precisely the Cu02 planes 

This first paper is devoted to the study of a single bipo- „ ^ ira • w i -i • -i r _l 

of cupratcs P|, it might exhibit similar features as more 

laron in the HH model in the adiabatic limit, assuming ... , , t , t • ^ 11 -,11 i • 

realistic models, in two-dimensional models with local in- 



classical phonons. Obviously the assumption that there 
are no quantum phonon fluctuations does not allow su- 
perfluid states (with many electrons). In the next paper 
jTzt , the quantum phonon correction to the adiabatic case 
will be studied. There, it will be shown that in some re- 



teractions, the bipolarons exist only for a large enough 
electron-phonon coupling and are always sharply localized 
(small bipolarons). We numerically calculate these bipo- 
larons by using a continuation method of these solutions 

from the anti-integrable limit |Q, where the electronic 

gions of the parameter space, there is indeed a drastic , r ■ , 1 • 

transfer integral is zero. 

reduction of the quantum bipolaron's effective mass due 

to quantum resonances between several almost degenerate The ground state of the bipolarons in this limit can be 

adiabatic bipolaron structures. A large part of the scien- easily found and consists of cither a bipolaron localized at 

tific material of these two papers can be already found (in a single site (SO) or of two uncoupled polarons at arbitrary 

French) in the PhD dissertation of one of us ||T8|. different sites, but there are many other states with larger 

energy that are combinations of singlet states (multisin- 

Some numerical studies of the bipolarons in the one- ,,nt,t ri 1-1 , 1 • 1 

glets). Many of these bipolaron states can be continued 

dimensional adiabatic HH model, were already presented , ; , , n ■ . 1 ■ r 1,1- 

when the transfer integral varies Irom zero and their en- 



ergies can be compared. Although the bipolaron (SO) or 
the singlet bipolaron (SI), persist with the lowest energy 
in large parts of the phase diagram, it is found that a 
quadrisinglet state (QS) becomes the ground-state in an 
intermediate regime of parameters. 



in ref. p9| (as well as few preliminary studies in two dimen- 
sions). Bipolarons always exist in one-dimensional mod- 
els as expected, but when the Hubbard term v increases 
from zero, a first order transition occurs between the single 
site bipolaron (SO) and a bipolaron (SI) composed of two 
bounded polarons on two neighboring sites in a magnetic 

singlet state. It was observed that the classical mobility of We show that we can reproduce quite accurately the 
the bipolaron (assuming the lattice dynamics is classical) same phase diagram by choosing variational wave func- 
was significantly enhanced in the vicinity of this transi- tions for the electrons made from simple combinations of 
tion. Owing to the presence of the Hubbard term, quite exponentials reproducing the main characteristic of the 
small bipolarons could become nevertheless highly mobile spin structure of the bipolaron. ( This is an extension of 
over hundreds of lattice spacings. the variational method used in ref. |^f| ) . Further exten- 



L. Proville, S. Aubry; Small Bipolarons in the 2-dimensional Holstein-Hubbard Model 



sions could be developed later for the many-body prob- 
lem. 

We investigate the properties of all the obtained so- 
lutions by calculating their binding energies, their pin- 
ning and breathing modes and also their Peierls-Nabarro 
energy barrier. We find a substantial softening of their 
pinning (and breathing) modes and a sharp depression 
of the PN energy barrier in the region where the (QS) 
bipolaron becomes the ground-state. Although the classi- 
cal mobility of the bipolarons never becomes as large as 
in the one-dimensional case p9[ |, it is sufficient to favor 
a good quantum mobility p7| in a specific region of the 
phase diagram. 

2 The Model 



dispersionless optical phonon branch with order of mag- 
nitude a tenth of an eV at most. 

g is the constant of the on-site electron-phonon cou- 
pling which may physically range from zero to a fraction of 
an cV. The on-site electron-electron interaction is repre- 
sented by a Hubbard term with positive coupling v which 
may range physically from negligible to large values of the 
order of 10 eV. 

Choosing Eq — ^g^ /Tiloq as the energy unit and intro- 
ducing the position and momentum operators: 



Pi = «T — [aj - a^) 

hujo 



we obtain the dimensionless Hamiltonian: 



(2) 
(3) 



To keep in mind the physical magnitude of the dimen- 
sionless parameters involved in our reduced model, let us 
first write the Holstein-Hubbard Hamiltonian with all its 
parameters measured in the original physical units: 



n=-T C+^C,,„+Y.^u:,{ata,) 

<i,j>,(J i 

+ Ygni{ai + a,) +Yvni^^n,^l (1) 



The electrons are represented by the standard fermion 
operators C^^ and Cj^a at site i with spin cr =| or [. Then 
T is the transfer integral of the electrons between nearest 
neighbor sites < i,j > of the lattice. In physical systems, 
its order of magnitude is usually measured in eV. 

and a; are standard creation and annihilation bo- 
son operators of phonons. htoo is the phonon energy of a 



2 ' 2 



<i,j>,a 



2 



(4) 



The parameters of the system are now: 



En 



T 



1 hujQ 

7 = li^) 



V 2g 



(5) 



The parameter 7 measures how "quantum" is the lat- 
tice. The BCS theory requires g << huiQ-. that is, large 7. 
We are interested in the opposite regime of strong electron- 
phonon coupling: that is, g larger than the phonon energy 
hujo- Then 7 becomes small. 

Thus the adiabatic approximation, which is simply ob- 
tained by taking 7 = 0, becomes valid in the strong elec- 
tron phonon regime. We shall assume this condition in 
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this first paper. Then {ui} commutes with the Hamilto- where A is the discrete Laplacian operator in the 2D lat- 

nian and can be taken as a scalar variable. For a given set tice (Z^)^ defined as = ^j-i '^j where j S {Z^)^ 

of {ui}, the adiabatic Hamiltonian are the nearest neighbors of i e (Z^)^. 

Unlike the singlet states, the eigenenergies of the triplet 

/ 1 I \ t states do not depend on the Hubbard term U since ipf^ = 

Had = Y^i l^uf + -u,n, + U n,^n,i j - - ^ C'+.Q, <t 

i ^ ^ <ij>,(T and thus are just the same as for noninteracting electrons. 

(6) 

commutes with the total spin of the system. 

Thus, the eigenstates of a system with two electrons 
are either nondegenerate singlet states or three-fold degen- 
erate triplet states. The wavefunction of the singlet state 
has the form 



Taking into account that in our model, the transfer in- 
tegrals with amplitude t > connect only the nearest 
neighbor sites, it is straightforward to check that the sin- 
glet state defined as ipij = I'^I.jl always has less energy 
than the triplet state with wave function {tpfA. As a re- 



suit, the ground-state of our system is necessarily a singlet 
<P>=Y^ C+T ^jU state with the form §) . 

The energy of (||) depends on {V'ij } and {ui} as 



where |0 > is the vacuum (no electrons in the system) and 
V'ij = "^j.i is normalized 



-^-<i>\A\^> (11) 



where the electronic density at site i is 



on the 2D lattice (Z^)^ (D = 2 being the lattice di- 
mension we consider in this paper). The wave function of 
the triplet state (oriented with the spin +1 in order to fix = ^^dV'jj l^ + iV'j.il^) (12) 



the ideas), has the form 



2J 



Extremalizing F{{'ipij}, {ui}) with respect to the nor- 
malized electronic state {V'ij} and the displacements {u^} 
yields the set of coupled equations (|l^) and 



where ipf, = —tpji is normalized and antisymmetric. 



Actually, the singlet wave and the triplet functions which n, ^ ~ q (13) 



are eigenstates of the adiabatic Hamiltonian both yield 
the same eigen-equation for their components tpij or tpf^ 



Fei{{ui}) is an eigenenergy of two interacting electrons 
in the potential generated by the lattice distortion {ui}. 
^Atfji^j + (^{"^i + ^j) + U6ij^ = Fel{{ui})^l;iJ Using eq.|l3|, the extrema of eq.|ll| are those of the varia- 

(10) tional energy 
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(14) 

for ipij normalized and where pi is given by eq.|l^. Then, 
it follows that 

- ^A%l;,^j + (y-\{Pt + Pj) + U5i^^ = Fei%l>i^j (15) 

and also that the solutions of this equation, the energy of 
the system is 

i^„({V,,,})=i^ei + ^E'^^ (16) 
i 

3 Numerical Continuations of Bipolarons 
from the Anti-lntegrable Limit 

3.1 Bipolarons in the Anti-lntegrable Limit 

In the anti-integrable limit t — Q, the adiabatic ground- 
state for two electrons is easily found. For U < 1/4, it 
consists of a pair of electrons localized at a single site i. 
This is the standard small bipolaron known in the litera- 
ture, denoted (SO) (see figj^)- For a bipolaron at site i, its 
electronic wave function is 

|iZ'>=C+C+|0> (17) 

and its energy is i^^ = ?7 — 1/2. 

When U > 1/4, the ground-state consists of two un- 
bound polarons localized at arbitrary different sites i and 
j and with arbitrary spins. It is thus degenerate and its 



the 2-dimensional Holstein-Hubbard Model 

energy F„ = — 1/4 is independent of the Hubbard interac- 
tion. When sites i and j are nearest neighbors, we define 
the bipolaron (SI) p5| , p^ (see figj^) with electronic wave 
function 

\^>-^(C+,C^^, + Cl,C+,)\9> (18) 

where i and j are nearest neighbor sites. 

Since a single polaron has the electronic spin 1/2, when 
the transfer integral t is small but not zero, a standard 
perturbation theory yields an antiferromagnetic exchange 
coupling 2t^ /U between the two spins of the uncoupled 
neighboring polarons. When the spins are chosen in the 
singlet state represented by eq.|l^, these two polarons have 
the energy F„ « —1/4 — /U. When they are not lo- 
cated at nearest-neighbor sites but at the lattice distance 
n, perturbation theory to order n yields an antiferromag- 
netic exchange coupling proportional to U{t/U)'^". Thus, 
for t << U, the minimum energy is obtained for nearest 
neighbor bipolarons in the singlet magnetic state (SI). It 
is maximum when U is close to and above 1/4, just when 
(SI) becomes of lower energy than (SO). For t fixed, it de- 
creases to zero when U increases. This binding energy also 
vanishes in the anti-integrable limit t. Unlike bipolaron 
(SO), bipolaron (SI) breaks the square lattice symmetry 
and is oriented either in the x direction or the y direction. 

When t is not very small, the spatial extension of the 
polarons goes significantly beyond single sites, and it is 
not obvious that a low-order perturbation theory holds. 
The true ground state might not be obtained by contin- 
uation of the solutions (SO) or (SI). There are infinitely 
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many other bipolaron states at i = (solutions of eq.p^, phonon coupling and only a small energy loss due to the 

which have been classified in appendix (^) |^ Some of Hubbard term (no doubly occupied sites). Moreover, the 

them are not very different in energy and could compete to peripheral electron can gain a substantial electronic ki- 

become the true bipolaron ground-state when t increases, netic energy by occupying N2 sites when A^2 > 1- In the 

Therefore, it becomes useful to test the ground-state of limit of N2 large, this energy gain can reach a maximum 



the bipolaron at t 7^ among the extrema of eq, 14 that at 2t. These bipolaron states have no continuous degener- 
are obtained by continuation from those calculated in the acy a.t t — and thus according to the implicit function 
anti-integrable limit at t = 0. theorem, they can be continued for t not too large. 

It is of course impossible to continue and to test nu- At t = 0, the electronic wave function of a star N2- 
merically the energy of all the solutions of eqs jl| at t = 0. singlet bipolaron centered at the origin is: 
The study of appendix (^) shows that the binding energy 

of the bipolarons is non-negligible only when the total y — ^ (^r!+ (j+ +(7+ (7+ )|0 > (19) 

number N ^ of occupied sites is not too large. The non- 

where ji, are neighboring sites to the origin. They could 

connected bipolaron states are discarded because at t = 

also be chosen farther away,but when the bipolaron be- 

they always have more energy than their connected com- 

comes too extended, its energy does not decrease suffi- 

ponent with the smallest energy, and at i 7^ 0, their ab- 

ciently to become the ground-state. We tested the most 

sence of connectivity is not favorable for gaining energy 

compact bipolarons which are star bisinglet bipolaron states 

from the electronic kinetic energy term with amplitude t. 

(BS) with N2 — 2 and ji, are the two neighboring sites to 
On the contrary, the star multisinglet bipolaron states ^^jgj^ direction x (or x and y), star trisinglet 

with one central site with electronic density pi = 1 (TVi - bipolaron states (TS) where 7V2=3 and 3, are three of the 
1) and 7V2 < 4 nearest neighbor sites with electronic den- neighboring sites of the origin, and the square symmet- 
sity p2 = I/N2 and energy F,, = -(1 + 1/N2)/S at t = ^-^ quadrisinglet (QS) (iVa = 4) which involves the four 
appear much more favorable for reducing their energy neighboring sites of the origin. 

when t increases. They are still spatially well-localized, ^^^g^^ -^^^^-^^ l^^^j^^^^ multisinglets with equal 

which allows an efficient energy gain from the electron electronic densities at the occupied sites might not be too 



1 Actually, this result should not be surprising since there high in energy and have been also tested. Although they 

are already infinitely many metastable states in addition to are continuously degenerate in the anti-integrable limit, 

the standard single bipolaron Q (see also section 5.4 in ref their degeneracy is raised when t ^ 0. We considered for 

[p^) in the pure Holstein model, which however never become example, the square symmetric quadrisinglet state (QS2) 

ground-state. which occupies the four corners ji, of an elementary square 
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of the lattice. One of its degenerate wave functions is where A/i is the normalization factor (chosen negative) 

— , 1 , , and K is some positive constant that we introduce to en- 

l'^>=E7g^>,T^tal0> (20) 

"5^"^' sure the convergence to a minimum energy state. Actually, 

with energy _F„ = —1/8 ., , , , , • ,, , • r 

it can be chosen to be zero m the domain of parameter we 

study. 

3.2 Numerical Technique of Continuation 

We find numerically that for n large enough, ifVi = 

The most efficient numerical techniques for the continu- T(!f'„_i) and its normahzation factor A/'„ converge to the 

ation of solutions of sets of equations as a function of a limits tf/ and TV, respectively, tf' is a solution of eq.|l^ with 

parameter, are usually based on a Newton method. For ex- the condition (|l^) and for the eigenenergy F^i = TV — 

ample, such techniques were developed efficiently for cal- K. This solution corresponds to the eigenvector of eqjl5| 

culating discrete breathers Q . In our case to calculate (where pi and pj are fixed) associated with the eigenvalue 

accurately adiabatic bipolarons on a 2D system, a rea- F^i which is such that F^i — K has the largest modulus. In 

sonable size should be 10 x 10. Then calculating the 10'' principle, the constant K is chosen large enough in order 

components of Vij with a Newton method, requires to that Fei is surely the lowest negative eigenvalue: that is, 

work with huge matrices containing 10* coefficients: that for the electronic ground-state. One can easily check in the 

is, to use a large-memory, fast computer. Actually, smaller anti-integrable limit that if = is an appropriate choice 

size conventional computers suffice if one uses appropri- when U < 1/2. Varying one of the model parameters by 

ate techniques needing a much smaller working space. This small steps, each solution is taken as a trial solution for 

technique does not allow to continue all solutions but only the next step. It is easy to determine whether the solution 

those which are locally stable (in particular, the bipolaron varies quasicontinuously or discontinuously. 

ground-state) and those that can be made stable by fixing solutions in the anti-integrable limit which are 

some spatial symmetries of the bipolaron. non-degenerate, it can be checked that the hypotheses of 

This method is quite simple in its principle. To solve implicit function theorem, are fulfilled. Thus continu- 

eq.[| with condition ^ , we start from a normalized trial ^^^^^^ principle possible []. For those which belong to 

solution of eq.^ = {0^^} with = cjij^i, and we ^ degenerate continuum, the conditions for applying the 

calculate recursively a new normalized trial solution tf'i = -.^^licit theorem are not fulfilled, but when some spatial 

T i^) = } as 

^ The implicit function theorem was already used in simi- 



lar anti-integrable limits, for example in ref. |24| for polarons 



J\fii/jij — A(f)i j + j U5i j — ^ ^ (0| + (j)^^ f.) ~ k \ (f>i j and bipolarons in the original Holstein model or in ref. p^] for 



(21) discrete breathers. 
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Fig. 1. Schemes of the bipolarons (SO), (SI) and (QS) appearing as possible ground-states 
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Fig. 2. Phase diagram of the bipolaron in the 2D Holstein-Hubbard model in the plane of parameters U and t. There are four 
phase domains separated by first order transition lines corresponding to bipolarons (SO), (SI), (QS) and two unbound extended 
electrons. Also shown are the limit of metastability of the bipolaron (SO) (dotted line), bipolaron (SI) (dot-dashed line) and 
bipolaron (QS) (dashed line). Insert: Magnification of the phase diagram around the triple point involving phases (SO), (SI) 
and (QS). 



symmetries or some constraints on the solution are fixed, 
the degeneracy at i = can lifted and this theorem ap- 
plies. 

In the anti-integrable limit, only (SO) for U < 1/2 and 
(SI) for < ?7 (and (Sn) with n > being the distance 



between two polarons) are numerically stable: that is, can 
be followed continuously from < = by using algorithm 
(^. Actually, we choose as initial solution at t = 0, the 
exact bipolaron solutions described above, which are (SO), 
(SI), (QS), (BS), (TS) and (QS2). Maintaining by force 
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the spatial symmetries of the solution at t = 0, the conver- laron. The region below this line is divided into three do- 

gence process becomes stable again, and the continuation mains separated by other first-order transition lines. For 

of these solutions is feasible. U small, the bipolaronic ground-state is (SO). When U in- 

The main advantage of our method is that it can be creases for t not too large, there is a transition line between 

performed on standard computers. Its flaw is that we might bipolarons (SO) and (SI). For larger i, this transition line 

not be able to follow continuously a solution that is mathe- bifurcates at a triple point at i sa 0.785 and U ~ 0.235 into 
matically continuable. Actually, rather few bipolaron states two first order transition lines which both join the transi- 

are continuable. In contrast, our method is very reliable for tion line with the extended states. In between the fork that 

finding the true bipolaron ground state, because it brings is generated, there is a small domain where the bipolaron 

spontaneously the bipolaron solution to a local minimum (QS) that was initially unstable for t small, recovers its 

of the variational energy. stability and even becomes the ground-state. Other bipo- 

Actually, we checked that when there is no symmetry laronic structures continued from the anti-intcgrable limit 

constrains and independant of the initial trial solution, at < = appear as minimum energy states in the domain 

in most cases our numerical algorithm converges sponta- shown on fig.^ The (QS) solution can be viewed as a lo- 

neously toward the same bipolaron state, which then can calized RVB state similar to that proposed by Anderson 

be considered as the true bipolaron ground-state. How- some years ago in the pure Hubbard model in 2D as 

ever, this situation does not occur in the vicinity of the a theory for superconductivity in cuprates. 

first order transition lines where we can obtain a few dif- ^ i i ■ /r~>o\ i • i i 

in our model, this (QS) bipolaron has the quantum 

ferent bipolaron states depending on the initial condition, j_ / \ ^ ii i • i- ^ • t i • 

° symmetry (s) because the kinetic energy term is Laplacian- 

but then their energies can be easily compared to find the ,.1 tt ^ j r j- /Fl\ ■ • ^ ii 

° J J- jjj^g^ However, the study of appendix (|y) m the anti-mtegrable 

ground state. limit, suggests that it is close in energy to other states with 

quantum symmetry (s') or (d). Such symmetries could be 

3.3 Bipolaron Phase Diagram „ , , ,. , ^ , , . ^. ^i r r , ■ 

° favored by slight model variations on the form of the ki- 

The ground-state for a pair of electrons is obtained by ^"^^i*^ energy. 

comparing the energies of many bipolarons continued At the triple point, the bipolaronic structure of our 



from the anti-integrable limit (see also |19 ). For larger model is degenerate between three states (SO), (SI) and 

i, the ground-state corresponds to a pair of electrons ex- (QS). FigjH shows the profiles of the electronic density 

tended over the whole system. There is a first order tran- for these three types of bipolaron, which have the same 

sition line, when t becomes smaller, at which the two elec- energy. Interestingly, they extend significantly over only a 

trons bind with each other and self-localize into a bipo- few sites, and thus can be called small bipolarons. 
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The binding energy of a bipolaron is defined as the dif- 
ference between the smallest energy state where the pair of 
electrons is unbound, and the bipolaron energy. Depend- 
ing on the parameters, this unbound state could be either 
two extended electrons with opposite spins in the plane 
wave state at zero momentum or two polarons localized 
far apart. The variation of this binding energy versus U 
and for several values of t is shown fig.§. 

At the triple point, the binding energy of the degen- Fig. 3. Profile of electronic density versus site i at the triple 
erate bipolarons (in that case, to produce two extended point t = 0.0779, U — 0.234 for bipolarons (SO), (QS) and (SI) 
electrons) is much smaller than the binding energy of bipo- along their symmetry y-axis and the transverse x-axis. These 
laron (SO), at the same value of t but at C/ = 0. However, three bipolarons have the same energy. 




the bipolarons (SO), (SI) and (QS). For that purpose, the 
variational forms have to be chosen appropriately under 



it still has a substantial value that is physically far from 
being negligible. 

The binding of bipolarons (SI) and (QS) are physically 

, , , . , i J 1 • r • ■ rni 1 • two conflicting constraints. On the one hand, they should 

better interpreted as being of magnetic origin. 1 hcsc bipo- 

1 1 • 1 i 1 1 1 11 ..Li be physically realistic enough in order to mimic the real 

larons can be viewed as two closely bound poiarons with i- ^ j o 

1 /r> mi • 1 • 1- • ,1 1 , ,1 • ground-state. On the other hand, the analytical calcula- 

spms 1/2. Iheir binding energy is mostly due to the spin 

energy gain obtained by lifting the spin degeneracy as a 

singlet state. 

In the vicinity of the triple point and specifically in 
that region, the quantum lattice fluctuations (7 7^ 0) will 
also lift the degeneracy between the three degenerate bipo- 
larons (SO), (SI) (in both directions x and y), (QS) re- 
sulting in a sharp mass reduction (or equivalently a large 
tunneling energy or a large band width) (see 0). 



tions of their variational energy should be practically fea- 
sible. 

In ref . |^ , it was shown that an exponential form cen- 
tered at the occupied site with a unique variational param- 
eter, was a good variational form for a single polaron, re- 
producing accurately its quantitative properties. We choose 
a similar normalized variational form for the electronic 
wave function of bipolaron (SO) located at the origin 



4 Variational Calculation of Bipolarons 



1- A2 



V;f° = aa(I*i+IjD with A = ( — ) 



(22) 



We now reproduce, with good accuracy, the phase diagram (for i — {ix,iy), we set \i\ — \ix\ + \iy\)- This variational 
shown in fig|| using simple variational approximations for form is easily extended to the electronic wave function of 
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The variational form for the electronic wave function 
of bipolaron (QS) centered at the origin is a combination 
of four of these variational forms in the four directions of 
the square lattice, but now it becomes useful to introduce 
two variational parameters A and /i instead of only one, to 
distinguish between the spatial extension of the polaron 
that is at the center from those that are the periphery: 
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Fig. 4. Binding energy versus U of bipolaron (SO) (thin dotted 
line), (SI) (thin dot-dashed line) and (QS) (thin dashed line) 
at t = 0.05 (top) (compared to two remote polarons) and t — 
0.08 (bottom) compared to two extended polarons. The upper 
envelope (thick line) is the binding energy of the ground-state. 
Insert: magnification at the first order transition between (SO) 
and (QS). 



bipolaron (SI) in a singlet magnetic state located at sites 
(0,0) and (1,0): 



''■^ \/2 



with B 



(l-A^) 



2\2 



(1 + A2)V1 + 6A2 + A4 



(23) 



QS 

c 



78 



(li.l + ljyl) 



± 



{\i^\ + \iy\) 



where for normalization 

1 + 



^(;^(bx±i|+bj,|) _|_ ;^(bx|+b«±i|)) (24) 



= ( 



[(1 + X^f + A2(3 - A2)(1 + A^) + 8A^ 



(25) 



(1 - Am)4 

The energy (|4|) can be analytically calculated with 
the variational forms (p2|), ( [2^ ) and (p4[). Extremalizing 
the resulting energy with respect to the parameters A and 
fj, yields the energies of bipolarons (SO), (SI) and (QS) 
with a very good accuracy. We do not reproduce here these 
tedious calculation. We also remark that this variational 
method allows one to compute the bipolaron structures 
even when they become unstable so that they cannot be 
numerically continued with our method. Comparing these 
variational energies allows one to produce a phase diagram 
that is very close to the exactly calculated one (see fig-H). 

However, it is worthwhile to mention that the varia- 
tional form ( p4[ ) of bipolaron (QS) may yield some arte- 
facts which are not found in the exact numerical calcula- 
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Fig. 5. Comparison between the exact phase diagram fig.^ 
(thick lines) and its approximate calculation (dashed lines) 
secj^ 

tions (as often occurs in variational calculations). Fortu- 
nately, they occur in parameter regions where this solution 
is not the ground-state, and thus do not affect the phase 
diagram. First, there is a first order transition of A and 
/i near the anti-integrable limit. Second, there is another 
anomaly when increasing U. It is found that ^jj^^ bifur- 
cates onto a unbound solution where ^ — I. This corre- 
sponds to an unbound pair of electrons in the spin singlet 
state where one electron is self-localized as a polaron and 
the second one is extended. As a result, the validity of the 
exponential form jg limited to the (QS) region, that 
is when this bipolaron is the ground-state. 

5 Internal Modes and Peierls Nabarro 
Barriers 

The phonon frequencies of the bipolaron can be easily cal- 
culated within the standard Born-Oppcnheimer approx- 
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imation (in units ^^7) as explained in p9[ |. It is found 
that the bipolarons exhibit several localized (or internal) 
modes. The breathing mode has the same symmetry as 
the bipolaron. The pinning modes are spatially antisym- 
metric and tend to move this bipolaron cither in the x 
direction or the y direction. Fig.^ shows the variations of 
their frequencies with U. It is found that in the region 
of the triple point where three bipolaronic structures are 
almost degenerate, both the breathing and the pinning 
modes, soften significantly (approximately by a factor 2). 

These weak frequencies for the internal modes can be 
considered as evidence that the self consistent potential in 
which the bipolaron is pinned becomes rather fiat, which 
means a small Peierls Nabarro barrier (PN) . It is thus use- 
ful to calculate precisely this PN energy barrier in order 
to confirm this conjecture. In addition, it is found that 
the paths that yield the lowest PN energy barrier vary in 
the parameter space. The several ways to move the bipo- 
laron are sometimes almost equivalent in energy. These 
paths should play a role in the quantum tunnelling of the 
bipolarons 

The PN energy barrier is the minimum energy that 
must be provided to the bipolaron to move it by one lat- 
tice spacing. For that we have to determine a continuous 
path of bipolaronic configurations which connects the ini- 
tial bipolaron to a shifted equivalent bipolaron. There is 
a maximum of energy along any path, and the minimum 
over all paths of this maximum (called minimax) yields 
the PN energy barrier. 
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Fig. 6. Phonon frequencies versus U of the pinning mode 
(thick line) and breathing mode (thin line) for bipolaron (SO) 
(dotted line), bipolaron (SI) (dot-dashed line), bipolaron (QS) 
(dashed line) at t = 0.05 (top) and t = 0.08 (bottom). When, 
these bipolarons are ground-state, the corresponding lines are 
plain. Vertical lines indicate the first order transitions. 

To move a bipolaron with electronic wave function 
{ipn ml from site i to a neighboring site j (where {-0^ — 

j^j-i^ra+j-i\)^ consider a continuum of bipolaronic 
solutions {'!/'n,m(c)} which depend on c for cq < c < ci, 
and such that {■(/'n,m(co)} = {■0,\,m} and {V'«,m(ci)} 

It is convenient for simplicity to choose as variable 
c(tf'), one of the bipolaron components or a simple func- 



tion of them. For any continuous path that connects the 
bipolaronic ground-state — {ipl^ ml at site i to the same 
configuration — {-0^ ^} at an equivalent neighboring 
site j, c must take all the values between cq = c{'F^) to 
c\ ~ c{'F^). For each value of c, the energy of the bipo- 
laronic state will be always larger than or equal to the 
minimum of energy of the bipolaronic configuration where 
the component corresponding to c(!f') is fixed to c. Thus, 
starting from the initial ground-state configuration, and 
following continuously this minimum by varying this con- 
straint c, we may pull continuously the bipolaron from one 
site i to its neighboring site j. 

For that purpose, the choice of c{'F) has to be appro- 
priate to obtain a path of bipolaronic configurations that 
connects continuously the two bipolaronic configurations 

and If"^ and that yields the lowest minimax. We guess 
intuitively that the bipolaron could be effectively pulled 
only if this constraint affects the " main body" of the bipo- 
laron instead of a minor component. For our investiga- 
tions, we found several continuous paths of configurations 
competing for providing the minimax. We obtain them by 
using several kinds of constraints for a bipolaron at site i 
moving to site j which may be: 



V'i.i = c 

i>i.j = -4^3,% = c 

0i,i - -^j,] = c 



(26) 
(27) 
(28) 
(29) 
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(fc 7^ i is a neighboring site of j and bond j — k could be 
either collinear with or orthogonal to bond i — j.). These 
constraints c are easily taken into account with few minor 
changes in the numerical programs described above min- 
imizing the variational form (|lj). When ipis, or ipj^k 
is fixed to c, it suffices to drop the corresponding equation 
(|l5|). When -04,1 — 'ipjj — c, it is convenient to define the 
new variable cj) — {ipi,i + V'jjO/V^- Then, the variational 
form (p^, depends on 0, c and ipn.m for (m, n) ^ (i, i) and 
7^ {jij)- The vector with components (j) and {V'm.n} (with 
(m, n) 7^ (i, i) and ^ has norm -^/f — which is 

fixed by the constraint c. Extrcmalizing ( p^ ) with respect 
to its free variables yields a set of equations that differ 
slightly from (^5|) , although they depend on c as a param- 
eter. They can be solved with the same iterative method 
as before. 

We may thus obtain a continuous path of configura- 
tions parameterized by c and connecting the bipolaron 
ground-state to the same state shifted by one lattice spac- 
ing. The extrema of the energy F„(c) given by ( |l4|) (which 
satisfy dF^ (c) /9c = 0) correspond to bipolaronic solutions 
without constraint. Actually, they can be identified among 
the bipolarons that are classified in appendix (^). When 
one of these bipolaron is found spatially symmetric with a 
symmetry center at the middle of the bond < i, j >, there 
is no need to continue the path beyond this point because 
it is clear that it can be completed by symmetry. We test 
the different constraints (^9|) and among those that yield 
a continuous path, the lowest maximum is considered as 
the minimax. 
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Fig. 7. Energy variation versus c of the bipolaronic ground- 
state state under the constraint -j/ji^i — -i/ij j = c. Bipolaron (SO) 
is initially at site i and j is the neighboring site towards which 
this bipolaron moves, t — 0.08 and (7 = (plain line), (7 = 0.1 
(long dashed line). Vertical lines indicate the location of the 
energy extrema corresponding to the initial bipolaron (SO). 
Insert: Variation of the Profile of electronic density along the 
continuous path of the bipolaron a.t t — 0.08, U = 0.1. 

The PN energy barrier is then the difference between 
the minimum of energy and the maximum along this best 
continuous path. 

When the PN energy barrier is smaller than or at most 
comparable to the binding energy of the bipolaron, the 
two polarons remain surely bound during their continuous 
lattice translation and it is then reasonable to believe that 
the minimax obtained with the above method is correct. 
However, there are regions in the parameter space where 
this condition is not fulfilled, and a precise determination 
of this PN energy barrier might be questionable. In any 
case, the obtained value, if not exact, is necessarily an 
overestimate. 
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When (SO) is the ground-state and U is not too large, ate unstable bipolarons (with one unstable mode), which 
the energy variation versus c = — J starting from the are nothing but the star sister bipolarons (SI/SO) de- 
ground-state bipolaron (SO) is shown in figj^. This path scribed in appendix ( [A.2| ) ^. Actually, this bifurcation line 
does not need to be continued beyond the point c = between (SI/SO) and (SI) appears on fig.^las the left bor- 
because it exhibits a minimax at c = corresponding to a der line of the domain of metastability of bipolaron (SI), 
spatially symmetric bipolaron and thus can be completed In that regime, the motion of the bipolaron involving the 
by symmetry. This bipolaron has only one unstable mode minimum energy consists in first stretching bipolaron (SO) 
and is the continuation at U small of the stable bipolaron into bipolaron (SI) along one lattice direction, and next in 
(SI) beyond its bifurcation point. At U — 0, its electronic squeezing this bipolaron in the same direction to recover 
state corresponds to two electrons with opposite spins in the bipolaron (SO) translated by one lattice spacing. This 
the lowest eigenstate of the potential generated by the lat- feature is identical to those found for the two-site model 
ticc distortion (Slater determinant). The analysis of ap- in appendix (|b|). 

pendix (g) suggests that for U < 0, this bipolaron can ^j^^^^ bipolaron (QS) (which does not exist for the 

be continued as the two-site bipolaron (2S0) in the anti- two-site model) becomes the ground-state instead of (SO), 

integrable limit. the PN energy barrier should be studied from this ini- 
tial configuration. Fig.0 shows the energy variation versus 



When U becomes larger, the previous constraint does 
not always yield a continuous path, and another constraint 
ipij = c is found efficient for providing a continuous path 
with a minimax. The energy variation versus c = Vij 
starting from the ground-state bipolaron (SO) is shown in 
figjl for some bipolarons. We observed that in addition 
to the minimum (SO), it exhibits two other extrema. The 
second minimum corresponds to the spatially symmetric 
bipolaron (SI). Again, there is no need to construct a com- 
plete path reaching (SO), since this path can be completed 
by symmetry. This figure shows that a pitchfork bifurca- 



Tpij = c starting from bipolaron (QS). The continuous 
path exhibits another minimum corresponding to the sta- 
ble bipolaron (SI), which is spatially symmetric. Again, 
the continuous path can be completed by symmetry. There 
is a minimax which correspond to another bipolaronic con- 
figuration, which wc did not analyze in detail but is likely 
to be the star trisinglet denoted (TS) described in ap- 



pendix (A.l). This curve also demonstrates that for this 
value of t, the bipolaronic ground-state changes by a first- 
order transition from (QS) to (SI) when increasing U. 



tion occurs for the minimax when U increases from zero » worthwhile to note a similar phenomenon observed 

(at fixed t). The unstable bipolaron (2S0) bifurcates into when narrow discrete breathers become mobile Interme- 

a minimum corresponding to the stable bipolaron (SI) diate discrete breathers breaking the lattice symmetry were 

and two symmetric minimax corresponding to intcrmedi- also found to appear |28| . 
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Fig. 8. Same as fig.^but with the constraint i/iij — c and dif- 
ferent parameters, top: t = 0.03, U = 0.1 (full line), U = 0.15 
(long dashed line), U = 0.2 (dashed line) insert: magnification 
bottom: t = 0.08, U = 0.18 (fuU line), U = 0.2 (dashed line). 
Vertical lines indicates the location of the energy extrema cor- 
responding to bipolaron (SO) and to bipolaron (SI). 
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It is also worthwhile to note that there is also a PN en- 
ergy barrier between the bipolaron (SO) and (QS), which 



have the same symmetry (see fig, 10). We tested that it 
does not generate any path with a lower PN energy bar- 
rier when shifting the bipolaron (SO) or (QS) by one lattice 
spacing. For that, we compare the energy barrier obtained 



0.2 0.4 0.6 

Fig. 9. Same as fig.^ but for bipolaron (QS) with the con- 
straint tpij = c at t — 0.08, U = 0.235 (dotted-dashed line), 
t = 0.075, U = 0.25 (dashed line) and t = 0.07, U = 0.26 
(full line), (Vertical lines indicates the location of the energy 
extrema corresponding to bipolaron (QS) and (SI). 

for shifting (SO) by one lattice spacing via the direct path 
(SI) Eb{SO SI), to that obtained by the indirect path 
(SO) via (QS) and (SI) involving the jump of two con- 
secutive barriers EeiSO QS) and Eb{QS SI). We 
found that the energy barrier between (SO) and (QS) was 
always relatively too high to favor the indirect path. 

When bipolaron (SI) becomes the ground-state, there 
are two PN energy barriers depending on the direction it is 
displaced, transversally or longitudinally. If it is displaced 
longitudinally in the direction of the bond (i, j) where (SI) 
is localized, the minimax may be obtained by varying the 
constraint ipj^k — "tpkj — c which tends to displace (SI) 
longitudinally. Fig.|ll| shows this energy variation versus 
c starting from bipolaron (SI). The maximum along this 
path corresponds to the longitudinal star bisinglet bipo- 
laron (BS). For t = 0.03 and U > 0.28, it costs less energy 
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Fig. 10. Top: Radial PN energy barrier versus tpi^i — c be- 
tween (SO) and (QS) at t = 0.08, U = 0.23 Insert: Profile 
variation of the bipolaron electronic density along the same 
path. Bottom: Same at t = 0.092, U = 0.1 between bipo- 
laron (SO) and the extended state at tpi^i — (full line) and 
t = 0.092,(7 = 0.204 between bipolaron (SO), (QS) as an inter- 
mediate state and the extended state (thin line). Insert: Bipo- 
laron profile (bottom left) and magnification of the (QS) region 
(top right). 

to use this path with minimax (BS') than to use the path 
with minimax (SI/SO) passing by bipolaron (SO). 

The transversal PN energy barrier of bipolaron (SI) 
can be also calculated. Actually, the transversal motion of 
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(SI) with the lowest PN energy barrier has to be done in 
two steps (in the anti-integrable limit). If we denote by 
(i,j, fc, Z) the corner sites of an elementary square of the 
2D lattice and move (SI) from the bond i — j to bond l — k, 
then (SI) rotates once by tt/2 around the center site i and 
again by 7r/2 but around the center site I. These two jumps 
have the same PN energy barrier. It can be measured from 
the height of the minimax determined by the constraint 
'4'i,i = = c. This path yields the 2 star multisinglet 
(BS') with a diagonal symmetry axis where "tpij — ipij and 
where the two branches and (i,l) are orthogonal. It 
is found that the longitudinal and transverse PN energy 
barrier are almost equal. 

The precise determination of the PN energy barrier be- 
comes more delicate close to the border line of the phase 
diagram figj^ with the domain of extended electrons. Here 
the binding energy of the bipolaron becomes very weak. 
To move a bipolaron, it may cost less energy to pass the 
energy barrier for breaking the pair of electrons (fig.^j), 
producing extended electronic states, and next to pass a 
new equivalent energy barrier for reconstructing the bipo- 
laron at another site. 

Figs.^ gather the resulting PN energy barrier versus 
U and for several values of t obtained by comparison of 
these different paths (for that reason we have broken lines 
with possible discontinuities). The essential result is that 
close to the region of the triple point between bipolaron 
(SO), (SI) and (QS), the Peierls-Nabarro energy barrier 
sharply drops and reaches the same order of magnitude as 
the binding energy of the bipolarons. When U slightly in- 
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Fig. 11. Same as fig.[7| but for bipolaron (SI) moving in the 
longitudinal direction with the constraint ipj^k ~ c or rotating 
transversally with the constraint tpi,i — c &t t = 0.01, U = 0.3 
(upper lines), t — 0.03, U = 0.3 (lower lines). (Vertical lines 
indicates the location of the energy extrema corresponding to 
bipolaron (SI) and (BS) or (BS')). Although almost the same, 
the curves for the transversal motion are slightly lower than 
those for the longitudinal motion. Insert; Variation of the Pro- 
file of electronic density along the axis i—j of bipolaron (SI) for 
the two continuous paths of the bipolaron at t = 0.03, u = 0.3 
corresponding to the longitudinal motion (left top) and the 
transversal motion (right bottom). 

creases beyond this point, the binding energy of the bipo- 
laron sharply decreases. Conversely, when U slightly de- 
creases, the PN energy barrier sharply increases. In that 
region, the paths aUowing a shift by one lattice spacing of 
a bipolaron with the smallest PN energy barrier involves 
the successive transformations . . . (QS) (SI) (QS) 

There are regions in the parameter space where the 
bipolaron binding energy becomes very small, for example 



t=0.03 
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Fig. 12. Peierls-Nabarro energy barrier (thick lines) and Bind- 
ing energy (thin lines) of ground-state bipolarons versus U for 
different values of t: t — 0.03 (top), t — 0.07 (middle) and 
t — 0.08 (bottom). The lines are broken because of the change 
of minimax. Inserts: magnification in the region of first order 
transition. 
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when t is small and U > 1/4) (e.g. fig, 12 at t = 0.03). of the bipolarons soften significantly in that region. More- 

The PN energy barrier of a bipolaron is then practically over, the Peierls-Nabarro energy barrier (FN) of the bipo- 

equal to the PN energy barrier of a free single polaron. In laron in that region is strongly depressed, which improves 

other regions, close to the first-order transition border line the classical mobility of this bipolaron. This effect is re- 

with the extended states, the bipolaron binding energy lated to the appearance of several intermediate metastable 

also becomes very small, but then the PN energy barrier is bipolaronic state which have almost the same energy. A 

practically that which has to be overcome for the electron small variation of parameters (t/, t) in that region suffices 

delocalization. either to lift the near degeneracy, with a PN energy bar- 
rier which grows very fast, or to depress sharply the bind- 

6 Concluding remarks ing energy of the bipolaron itself. The energy landscape 

around the bipolaron has been explored. It has been found 

The interplay between the electron-phonon coupling and 

that it is quite flat in the region of the triple point with 

the direct electronic repulsion has been treated accurately 

several minimum energy states close in energy and small 

in the adiabatic Holstein-Hubbard model in two dimen- 

energy barriers between them. 

sions. Numerical investigations complemented by analytic 

... Ill,- • ij 1 T r These features strongly support the coniecture that 

variational calculations yield the phase diagram of the j ft- j 

- ii-i 1-1 -ii- the quantum tunneling of the bipolaron will be strongly 

ground-state of a single bipolaron, which consists of sev- ^ o t a j 

, , . i 1 1 i2 i 1 i -i- T / enhanced in the vicinity of this triple point due both to the 

eral domains separated by nrst order transition lines (see f f 

c ^^ Ti • r J i-1 i- j-£r iu- i • j- ^ i-i ^ Small PN energy barrier and to the hybridization between 
ngJ3) . It IS found that the diflerent bipolaronic states that 

, , . J , , ■ , , , 1 i- ■ i 1 1 1- -i J the nearly degenerate states. This assertion will be con- 
are obtained, already exist at the anti-mtegrabie limit and ° 

, i ir ii • 1- -i 1 i- i- firmed by the results of the next paper where the quantum 
can be generated irom this limit by continuation. 

mi . • i X- • • ii 1 1- lattice fluctuations will be treated as perturbation through 

ihere is an interesting region m the phase diagram ° 

. jiiu a tight binding model mTl. 

where the bipolaronic ground-state becomes a quadrism- ° ° i^ 

glet bipolaron, which is a superposition of four singlets Unlike the conclusion of ref . |TT| , we find a plausible 
sharing one central site. The binding energy of that bipo- mechanism for a drastic reduction, under specific condi- 
laron is the result of the spin resonance between a strongly tions by several orders of magnitude, of the effective mass 
localized polaron and a peripheral electron localized on its of a bipolaron while preserving a relatively large binding 
nearest neighbors. energy. Let us recall that fig.^ shows that the binding en- 
There is a triple point where the three kinds of bipo- ergy close to the triple point is still about O.OOSE'q. Since 
laron coexist with the same binding energy, which is still Eq — 8g'^ /hujQ could reach in some realistic physical mod- 
significantly large and non-negligible. The internal modes els a magnitude of about lOeV, this binding energy can 
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be close to O.OSeT^ which corresponds in temperature units symmetry, has necessarily the trivial quantum symmetry 

to 500 K!. In the same region, the Peierls-Nabarro energy (s). However, it is not physically unrealistic to assume 

barrier has a value almost equal to this binding energy, that the electronic kinetic energy terms in Hamiltonian (|^) 

It is drastically reduced compared to what could be ex- might be different from a discrete Laplacian form ^. When 

pected for small bipolarons in standard theories. Note that there are second-nearest-neighbor electronic transfers with 

is about 50 times smaller than the Peierls-Nabarro energy significant amplitudes (but not necessarily as large as the 

barrier of the bipolaron (SO) a,t U — and the same value nearest-neighbor integral) and with appropriate signs, the 

of t. ground-state electronic wave function {ipi.j} of bipolaron 

j-u i i- f 4-1 u 1 iu- (QS) should have a (d) symmetry (see appendix O. A 

When the temperature of the system goes below this ' \ / j j \ ft- lj/ 

, c J superfluid state of such bipolarons with a degenerate in- 

characteristic temperature where bipolarons can form, and c- o 

^, , ,,. , , 1 111 ternal quantum symmetry could be perhaps related with 

II the tunnelling energy could reach comparable values ^ ^ 

/ , . , .,, , , . , , , N ii • rr i • 1 i the now well accepted fact that the superconducting order 

(which will be shown m the next paper), this ettect might r- o 

, rc ■ , , r a-jiii rvT^ • J- ■ parameter of cuprates has a (d) wave symmetry |29|l . Fur- 
be surfacient to favor a supernuid state at OK against ei- \ / j j il^sj 

, . , , . i- 1 1 ther works will investigate consequences of this quantum 

ther a bipolaron ordering or a magnetic polaron order- o i i 

rrn • i i 11 ■ i i 11 1 i symmetry. We expect that when the (d) wave symmetry 

mg. Ihis state could persist to unusually large temper- j j \ / j j 

mi ■ r ii 1-i- 1 • 1 • is favored by appropriate terms, the stability domain of 

atures. Ihere is, of course, another condition, which is j ± j- ; j 

j^ij^^iT ,1-1 ^. , the bipolarons that can take advantage of this symmetry 

that the direct bipolaron interactions are not too strong. a j j 

rpi. . 1-1 1 i u i. 4- u i-cciT u ii will be extended: that is, those of bipolaron (QS). 

I his IS very unlikely to be true at hall-faUmg, where the j 

polarons are close-packed, but this condition might be- principle, the method used in this paper for cal- 
come fulfilled when the density of electrons moves suf- culating adiabatic bipolarons could be extended to more 
ficiently far from half-filling. These quantitative results complex and realistic models. There are many kinds of 
in the Holstein-Hubbard model yield a more quantita- bipolarons in the anti-integrable limit, as shown in ap- 
tive support to earlier but less specific conjectures that pendix (g). It is not obvious that only bipolarons (SO), 
high Tc superconductivity could be explained by a well (si) or (QS) are competing as ground-states. More gen- 
balanced competition between electron-electron repulsion erally, if there are more transfer integrals between further 
and electron-phonon interaction neighbors, other N2 star multisinglets with iVa > 4 (eg. 



The methods used above should also work with other 



perturbations from the anti-integrable hmit. In the present * A reduced Holstein-Hubbard Hamiltonian on the copper 
paper, the Laplacian form for the kinetic energy implies square 2D sublattice of cuprates should involve more than 
that the bipolaron ground-state when it has the square nearest-neighbor transfer integral due to the oxygen bridges. 
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N2 = 8 etc. . . ) might become more favorable as bipolaron bridge for its hospitality during the completion of this 

ground-state in some cases. manuscript. 

Of course, the present study with only two electrons 

is by far not sufficient to describe real cuprates where the _ 
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other solutions with the same energy, which are simply ob- A.l Connected Bipolaron States with no Doubly 

tained by arbitrary permutations of the sites of the lattice Occupied Site 

If "01,1 = for any i, the electronic wave function is a 

For each solution, the set of occupied sites ieS'is t, ,• • r . •; - i;;; t n ^ 

normalized combination ot two sites singlet states dehned 

defined by the condition p; ^ 0. We call link a pair of . r-ri , , , , , , . 

m eq.|18|, and consequently these states and their energies 

sites (i, j) such that ibj i ^ ( which implies i €z S and , , , 

^ ^^,1 T- \ 1 do not depend on t/. 

i G S). A bipolaron state at i = is said to be connected at • i n i .i ; ; i i r 

Ns = JMi + N2 IS denned as the total number ot occu- 

if the graph generated by all the links is connected. • i •; at • ,i 1 r •; ■ n -jt ^ 

pied sites. A*! is the number ot sites m bi with density pi 

We first investigate the connected states of eqj30| which and the number of sites in S2 with density p2- Since 

imphes that when Vij 7^ 0, the total number of electrons is two, 



Nipi + N2P2 = 2. (32) 

Ud,,j^]ip, + p,)^F,i. (31) 

^ It follows from eqjlj that 

Fei is independent of the pair of connected sites Fv = ~g ^^P? ~ ~g(-^iPi + ^2^2) (^3) 

i 

Considering two different sites i and j connected to a third 



site n, it comes out that pi + Pn ^ Pj + Pn, which implies 

Pi — Pj. More generally, two occupied sites connected by -^v — -^\pi P2 ) — 



and equivalently from eq.^ and eq.pl| 

1, si 



= -7(Pi + P2) + ^{Nipl + N2PI) (34) 



some path with an even number of links have necessarily tj • ^ ^CT^ j /fT71\ j 

^ identilymg the two results ( |33| ) and (g4p, and using 

the same electronic density. As a result, the set of occupied rjT] ^ , ^. ^ i - 1 n 1^ 

■' ' ^ eq]3^, two solutions come out which are nrst 

sites S is the union of two disjoint sets of sites and S2 
where the electronic densities are the same. For i e S'l , 

1 1 
Pi — -rr and p2 — -rr with (35) 
the electronic density is pi = pi and for j Q S2 , Pj = P2 • ^1 ^2 

Moreover, sites i G Si are only linked to sites j G 52 and ~ ~ SN1N2 ^"^^^ 

and second 



vice versa. 



We consider separately the connected bipolaron states 2 

Pi = P2 = -rr with (37) 

J V c 



without and with doubly occupied sites i. These sites are 
defined by the condition ipi^i ^ 0. 



F. ^ (38) 



In the first case(|35|), we have pi 7^ p2 when A^i ^ N2. 
Then ipij 7^ when i G Si and j G ^2 or vice versa. 
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This condition determines a rectangular A^i x N2 matrix. 
The square of its real positive coefficients fulfills the linear 
eqs.^ 

E i^h - ^ (40) 

A particular solution of this set of equations is ipf j — 
l/{2NiN2). However, there are A^i + N2 linear equations 
to determine the N1N2 coefficients. Except for the case 
iVi = f or equivalently N2 — 1, (we have A^i ^ N2), 
this solution ipij is degenerate and belongs to a nonvoid 
bounded and compact domain defined by the positivity of 

The solutions with iVi = f appear especially interest- 
ing, not only because they are not continuously degenerate 
but because their energy Fy = -(iVs + l)/(8iV2) < -f/8 
is significantly lower than zero. It is not far above those of 
the bipolaron (SO) which is F„ = — 1/2 + [/ and those of 
the singlet bipolaron (SI) which is F„ = — 1/4 (see fig.p^ . 
We call them star multisinglets. (SI) and (QS) are star 
multisinglcts with N2 = I and N2 = 4: (see fig.|l|). 

In the second case (|3^) or in the first case when A^i — 
N2, the electronic densities at the Ng = 2Ni occupied 
sites are equal. Then there are in general Ns{Ns — f)/2 
nonzero coefficients V'i.j = "^j.i since there are no dou- 
bly occupied sites (V'i,i = 0). They fulfill Ng equations 
'^j'^tj ^ V(2-^s) (p^. Again, this system has a trivial 
solution, which is ■(/'f j = l/{2Ns{Ns - I)) for z 7^ j e S". 
However, when Ns > 4, this system of equations becomes 
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Fig. 13. At t = 0, energy as a function of A^'^: (a) star sisters 
solutions for U — 0.25, (b) star sisters solutions U = 0.4, and 
(c) star multisinglets (SI) for Ns = 2, bisinglet (BS) for Ns = 
3, trisinglet (TS) for Ns = 4, quadrisinglet (QS) for Ns = 5. 

underdetermined and yields continuously degenerate so- 
lutions that belong to a compact domain since again ipf ^ 
must be found positive. 

It is worthwhile to remark that although these so- 
lutions were assumed to be connected, when they form 
a continuously degenerate set this set may contain non- 
connected states just at the border of the compact do- 
main of solutions. For example, in this second case, let 
us split the set S of occupied sites in two disjoint sub- 
sets Ti and T2 with Mg > 2 and Ns - Ms > 2 sites 
respectively. Let us set ■(/'fc,; = fpi,k — for k £ Ti and 
I E T2, which corresponds to Ms{Ns — Ms) conditions. 
When Ns{Ns - I)/2 - Ms{Ns - Ms) > Ns, the equation 
Eji'h has the solution = f/V(Ms - l)Ns 

for i ^ j e Ti and ■(/',: j = ^/Vi^s -Ms- l)Ns for 
i 7^ J G T2. This situation is found to occur for Ns > 6. 
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A. 2 Connected Bipolaron States with Doubly According to eqjl^, the set of electronic states ipij — 

Occupied Sites at U ^ satisfy 

Let us require iV^ > 1 to avoid the onsite bipolaron (SO). El. — -|_ ^2 j e 5^ (47) 

Such connected solutions {ipi j} of eq.^ have at least one p„ 

y = "("h fo"^ 3 ^ ^2 (48) 

doubly-occupied site k ( i.e such that ipk,k 7^ 0). The set 

of occupied sites can be split in two sets Si and S2 with ^here pi and p2 are given by eqs.||. There are N1+N2 lin- 

electronic density pi and p2, respectively Let us consider ear equations for calculating iVi(iV2 + 1) positive numbers 
a doubly occupied site k which we assume to belong to ^it^^ « ^ '^'i and ^pf ^ = Vf,, with i e Si and j e 82- A 

the set of sites Si with electronic density pi. Let us also particular solution is obtained when all i'l^ = P2/{2Ni) = 

consider a site I € S2 such that Vfc,; 7^ 0. There exists (l-2C/iVi)/(iV,7Vi) and aU V^^, = {Nipi~N2P2)/{2Ni) = 

such a site since the solution is connected. Then, because (^1 " ^2 + 4f77ViA^2)/(A^iA^2). The positivity of Vf,, re- 

of eq.|l|, we have quires 

1111 -^TT^ < U 49 

Fei = --Pk + U = --{pk+Pi)^'-{pi+P2)=-^Pi + U ^NiN2 

(41) Actually, there are no positive solutions at all to eqsj47| 

which implies and ^ when this condition is not fulfilled. 

jj ^ Pi - P2 When conditions and (|9|) are fulfilled and when 

.,, ,, 1 , , . , , , , , ri mi A^i > 1, the number of variables exceeds the number of 

All the doubly occupied sites must belong to Ai. Ihe 

1 1 1 Ti 1 1 • 1 • T 1- i- i equations, and there is a compact set of degenerate posi- 
nonzero Hubbard amplitude U obviously implies distinct 

, , . , r 111 -1 1 111 five solutions to the linear eqs.E7^ and n§i. 

electronic density tor doubly-occupied and non-doubly- ' — ' ' — ' 

When iVi = 1 and when 



occupied sites. It follows from eq.|3^ that 

I + 2UN2 , l~2UNi ^2-1 .,^1 

PI ^2^^ and p2^2^^ (43) < U < - (50) 

which implies there is a unique solution to this set of linear equations. 

~ '2N2 ^ ^ ^ 2iVi ^'^^ ^ given Ns, it has the lowest energy Fy (plotted on 

in order that both pi and p2 be positive. fig-lll for different values of U). This solution is called 

We now obtain from eqs.p^,p^at t = N2 star sister. It can be interpreted as a mixing between 



1,9 9,1, , , , the bipolaron (SO) and the N2 star multisinglet. We de- 

Fv = ^{Nipl+N2pl)--ipi+P2) (45) i V ; 

note (SO/SI) the one star sister that mixes both (SO) 
and (SI) etc. . . . According to the implicit function the- 



4 

and substituting ^ in ^ 

1 ^■^r2 



IN AfiN2 + 2U{Ni N2) 1] (46) orem, this nondegenerate solution can be continued to t 
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nonzero except at the bifurcation points at U =^ 1/2 and 

U = (A^2 - 1)/4A^2 where this solution bifurcates with |;^^^|2 = iFei (53) 

(SO) and the A''2 star multisinglet, respectively. In the anti- 

If component a has doubly occupied sites (see ap- 

integrable limit t — its energy is larger than both those 



pendix. |A.2[) , then eq.^ implies 
of (SO) and A^2 multisinglet (see part |^ for t > 0). 

I A„ p = -NfF,i - U{N^ - TVf ) (54) 

A. 3 Non Connected Bipolaron States 

There are constraints for solving the second equation 
Let us now assume that we have a solution j } of eq.|^, because of inequalities (||) and , which imply 
which is not connected. Then, it can be decomposed into 

a sum of normalized connected components {V'fj}- We 1/7 1 

~ 2N°' \X p ^ 2iV" ^'^^'^ 
define Sa as a set of lattice sites that are connected with ^ ^ 

< (56) 

each other by a sequence of links, tpfj is proportional 4A^f AT^ l^ap 

to J for I e Sa and j G 5^ and zero elsewhere. The Conversely, we can construct non-connected bipolaron 
proportionality coefficient A„ is chosen in order that { V,^, } states which are combination of non-overlapping connected 
normahzed. Then we have states. Then the amplitude |A„|2 is defined by eqs. || or 

but then we have to choose F^i in order that the nor- 
ijj^ J — Xgipfj (51) malization condition ( |5^ ) is fulfilled. This is easy to do 

a 

when only components a with no doubly occupied sites 

with 

are involved. Otherwise, we have to take into account the 
y|A„p = l (52) „ „ 

^-^ constraints (M) and (Hy). There are many such solutions 



Two components a and (3 have no common occupied site, but we did not investigate them in detail. Some of them 

Then, {^l^} is a connected solution of eq.|§, for the ^re easy to find: for example, the nonconnected solution 

eigenenergy F^^ = F,i/\K\'' and the Hubbard term C/„ - "^^^^ components involving the bipolaron (SO) located 

Ij |2 on two adjacent sites. This solution is called (2S0) and has 

If component a has no doubly occupied sites, it does energy U 1/4. 

not depend on U. Then if TVf + = represents the '^^'^ ^e checked that the energy of the disconnected 

number of sites in each of the two groups of sites with ^^^^^ ^^^^^^ l^^S^^ ^h^" t^"^^ its components with 

different electronic densities defined at the beginning of smallest energy, 
this appendix, eq.pl] implies 
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B Two-site model 

It is instructive to analyze all the extrema of the varia- 
tional form (|lj) on a lattice reduced to only two sites i 
and J, because it can be explicitly calculated in all detail. 
However, a limitation of this restricted model is that in 
addition to the absence of extended states, the bipolaron 
(QS) cannot occur with only two sites. 

Setting ipij — X, ipj,j — U and i/^ij = z, using the 
normalization + + = 1, (|lj) becomes 

^V2t{x + y) Vl - x^ - y2 (57) 

Considering as equivalent the extrema obtained by sym- 
metries {x —t —X, y — + — y) and (x — )■ y, y — > a;), there are 
up to 4 kinds of extrema to this variational form at t = 0: 

— Bipolaron (SO) with energy U — 1/2 is a local minimum 
for U < 1/2 and becomes a saddle point with only 
one unstable direction for U > 1/2. There are two 
symmetric such solutions located either at site i or at 
site j. 

— Bipolaron (SI) with energy —1/4 is a local minimum 
for U > and becomes a maximum (with two unstable 
directions) for U < 0. 

— Bipolaron (2S0) is a non connected state consisting 
of ipi^i = tpjj ~ l/\/2 and ^Ji.j = i^j.i — 0. It is a 
maximum (two unstable directions) for U > and a 
saddle point (one unstable direction) for J7 < 0. 

— When < ?7 < 1/2, there is another extremum which 
is the 1 star sister (SI/SO) described in appendix 
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It is a saddle point with one unstable direction. It bi- 
furcates with bipolaron (SI) at J7 = and with bipo- 
laron (SO) a±U — 1/2. There are two symmetric such 
solutions located at site i or j. 

The minimax corresponding to the PN energy barrier 
for moving either bipolaron (SO) or (SI) is nothing but 
the unique saddle point which could be (2S0) {U < 0), 
(SI/SO) {0<U < 1/2) or (SI) (1/2 < U). 

At t = and [/ = 0, bipolarons (SI) and (2S0), which 
are both spatially symmetric, have also the same energy 
and the same electronic density (see fig.^). When t 7^ 0, 
this degeneracy is raised as shown in fig.p^. 

C Quantum Symmetries of Bipolarons 

Eventhough no nontrivial quantum symmetry appears for 
the bipolaronic ground-states of our model, it is worth- 
while to going forward now some further work of ours and 
discuss the possibility of nontrivial quantum symmetries. 
Actually, such symmetries are already latent in the anti- 
integrable limit and could be manifested easily in appro- 
priately modified models. 

As we pointed out, in the anti-integrable limit, only the 
modulus of ^prn,n IS determined but not the phases. This 
degeneracy is expected to be lifted by the perturbation 
from this limit due the electronic kinetic energy. However, 
it might not be completely lifted in some cases. 

This situation may occur for bipolarons associated with 
a lattice distortion (or equivalently an electronic density) 
which has the square symmetry of the lattice (group C^y). 
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Fig. 14. Bipolaron energies versus U at the anti-integrable 
limit t = (top) and at t = 0.05 (bottom). Ground-state (thick 
full line), stable bipolaron (full line), minimax (thick dashed 
line), maximum (thin dashed lines) 



This symmetry group has only two generators, which are 
for example the 7r/2 rotation and the reflection with re- 
spect to the X axis. Any of the symmetry transforma- 
tions change only the phase of the electronic wave function 

{^m,n} but not its modulus |'(/'m,ri|- 

There are three possible group representations for C4J, 
usually denoted A, B and E in textbooks of crystallogra- 



phy 1 31 . We denote them (s), (s') and (d) respectively^, 
(see fig.|l|) 

— When {i/'m.n} has the (s) symmetry, it is unchanged 
by any symmetry operation. This representation has 
obviously dimension 1. 

— When {^,n,„} has the (s') symmetry, {V'm,n} is changed 
into {— V'm.n}, by a 7r/2 rotation of the lattice. It is un- 
changed by reflection with respect to the x axis. The 
other transformations are obtained by combinations of 
these ones. This representation has also dimension 1. 

— For the (d) symmetry, {V'm.n} is changed into {z'i/'m.n}, 
by a 7r/2 rotation of the lattice and {t\}m.n \ is changed 
into l?/'^ „}, by reflection with respect to the x axis. 
This representation has dimension 2. 

We now note that, in the anti-integrable limit, bipo- 
laron (SO) always has the symmetry (s). For bipolaron 
(SI), which does not have the square symmetry but only 
an axis of symmetry, the symmetry is too poor to generate 
a d symmetry. Bipolaron (QS) is more interesting because 
it has the square symmetry for its electronic density but 
its quantum wave function may have three different quan- 
tum symmetries (s), (s') and (c?) respectively (see figjisl). 
These three states are degenerate in the anti-integrable 
limit but the electronic kinetic energy lifts this degener- 



In principle, the (d) symmetry characterizes the dimension 
5 representation I = 2oi the continuous rotation group 0(3) in 
three dimensions. We still use this terminology for the symme- 
try group of the square lattice although it is not appropriate, 
but because it is the most standard one in symmetry theory of 
superconductivity. 
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have nontrivial quantum symmetries (s') or (d) instead of 

(»)- 



1'0n= -l/VS 
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Fig. 15. Three (QS) bipolarons at the anti-integrable limit 
with the different quantum symmetries (s) (top left),(s') (top 
right) and (d) (bottom) . 



acy. For a Laplacian-like kinetic energy as in the model 
treated where < > in the prefered quantum symme- 
try is (s). However, symmetry (d) can be easily favored 
when there are next nearest-neighbor transfer integrals 
with negative signs in the electronic kinetic energy. 

In summary, we described in appendices (^) and (|c|) 
a systematic method that allows one to construct all the 
possible bipolaron solutions existing at the anti-integrable 
limit t = 0. There are bipolaron states with no continu- 
ous degeneracy and others with a continuous degeneracy. 
Bipolarons (SO) and star multisinglets with N2 small ap- 
pear to be the best candidates for bipolaronic ground- 
states when the electronic kinetic energy is switched to be 
non-zero. Star sisters may appear as minimaxes but are 
not found as ground-states. There are bipolarons with the 
square lattice symmetry (for example (QS)) which may 



